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Abstract 

A new criterion for A-stability of peer two-step methods is presented 
which is verifiable exactly in exact arithmetic by checking semi-definite- 
ness of a certain test matrix. It depends on the existence of two positive 
definite weight matrices for a given method. Although the initial ap¬ 
proach is different using properties of the numerical radius the criterion 
itself resembles the one from algebraic stability of General Linear Meth¬ 
ods. Known numerical algorithms for the computation of the unknown 
weight matrices suffer from rank deficiencies of the test matrix. For s- 
stage peer methods of order s — 1 this rank defect is identified with an 
explicit block diagonal decomposition of the test matrix in trivial and def¬ 
inite blocks. In the design of methods its coefficients are unknown and 
an explicit parametrization of A-stable peer methods of order s — 1 is 
presented with a weight matrix as parameter. This leads to a general 
existence result for any number of stages. The restrictions for efficient 
L-stable peer methods like diagonally-implicit and parallel ones are also 
discussed and such methods with 3 and 4 stages are constructed. 

Key words. Peer two-step methods, A-stability, numerical radius, general 
linear methods. 


1 Introduction 

A-stability requires that numerical solutions for the Dahlquist test equation 

y'{t) = Xy{t), A € C_ := {z € C : 5RA < 0}, (1.1) 

be uniformly bounded for t > 0. Although this concept seems rather limited 
addressing only the simple scalar equation it gains greater significance by the 
theorem of von Neumann [ig, 0 Th.IV.11.2] which is applicable to A-stable 
one-step discretizations of linear systems y'{t) = Ay{t) if A has non-positive 
logarithmic norm, e.g. 0. For Runge-Kutta methods the stability function R : 
C I—>■ C is a scalar function. Establishing A-stability here is possible through the 
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maximum principle by considering the E-polynomial \Q{ir])\'^ — \P(ir/)\'^,r] G K, 
of R{z) = P{z)/Q(z) on the imaginary axis if R has no poles in C_, [8]. This is 
an analytic criterion which may be difficult to check exactly for higher degrees. 
However, Hairer, Tiirke and Scherer also derived a purely algebraic criterion 
for A-stability of Runge-Kutta methods depending on the semi-definiteness of 
a certain matrix, H, m- It is a weaker version of the criterion for algebraic 
stability since the diagonal matrix of quadrature weights may be replaced by an 
unknown semi-definite weight matrix. For s-stage peer two-step methods the 
stability function is a matrix function M : C In this case we define 

A-stability (i.e. A-acceptability of the stability matrix) by 

e(M(z))<I, AgC_, (1.2) 

where g denotes the spectral radius. A matrix valued version of the von- 
Neumann theorem is due to Nevanlinna, m- It uses assumptions on the spec¬ 
tral norm ||M( 2;)||2 not on the spectral radius. However, it will be shown that 
such bounds follow from the new algebraic criterion. 

Peer two-step methods were introduced in |20j in a linearly-implicit version. 
They are time stepping schemes using s stages Ymi G K”, i = l,...,s, on 
time intervals with general stepsizes = tm+i — tm- However, the 

discussion of A-stability is restricted to a fixed stepsize = h, here. The stages 
Ymi are associated with the ODE solution at off-step points tmi = tm + hci with 
fixed nodes c^, i = 1,..., s. For an autonomous problem y'{t) = /(?/(t)) with 
right-hand side / : K" —>• R" and function evaluations Fmi = f{Ymi) the peer 
two-step method is given by 


Yfni 'y ^ 9ijFmj — ^ ^ ^ijYjjri—i^j hjyi ^ ^ dijPm-l^j-, ^ — 1, . . . , S. (1-3) 

i=l i=l 1=1 

The elements bij,gij, i,j = 1,..., s are the coefficients of the method and 
will be assembled in square matrices A = (a^), B = {bij), G := {gij) G 
For a constant stepsize these matrices are independent of the time step. For 
general dimension n G N stage vectors and function evaluations may be collected 
in matrices Y^ = (Ymi ,..., Fms) £ R®^”, = (/(^mi))yielding a compact 

representation of the method (11.31) as 


Yjn h^riG^F^ — BjyiY^—i -(- hyjiAjYiFjYi—\. 


(1.4) 


The structure of the matrix Gm determines the numerical effort of the scheme. 
In order to keep the effort of implicit methods moderate diagonally-implicit 
methods with a lower triangular matrix G may be considered [2] or parallel 
methods where G is diagonal, [501 HH US]- For nonlinear ODEs (11.41) is an 
implicit scheme. The notion ’peer’ refers to the fact that all stages Ymi have 
the same order. Hence, a polynomial predictor is available to obtain a linearly- 
implicit method of full accuracy by performing only one Newton step, |23j . 
Finally we note that peer methods are a special subclass of General Linear 
Methods (GLMs) where all internal stages are passed to the next interval. 
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Verification of A-stability for GLMs is more difficult than for Runge-Kutta 
and multistep methods. Floating point computations of eigenvalues are not fully 
reliable since g(M(0)) = 1 by preconsistency (see below). The Schur criterion 
for polynomials may be used, see mm. requiring global bounds on certain 
polynomials on C_. Verification may also be difficult here for higher degrees. 
In this paper a sufficient criterion is presented where only semi-definiteness of a 
certain (2s) x (2s)-matrix M (called a ’test matrix’) has to be checked. This test 
can be performed exactly e.g. by using rational arithmetic and a Cholesky-type 
decomposition if the matrix elements are rational. The criterion is based on a 
reformulation of the eigenvalue problem for M( 2 ;) and uses semi-definiteness in 
an inner product described by an unknown symmetric positive definite matrix 
Z (called a weight matrix). Using well-known properties of the numerical ra¬ 
dius of matrices an equivalent characterization by semi-definiteness of a certain 
block matrix M depending on two unknown matrices Z and W is obtained. This 
criterion, which is a feasibility problem of semidefinite programming f[25jl. is 
related to a reformulation of algebraic stability presented in [TO] . As in the case 
of Runge-Kutta methods the restrictions on the unknown weight matrices Z, W 
are weaker for A-stability. The construction of such feasible matrices for given 
algebraically stable GLMs has been discussed by Hewitt and Hill lilTO] and 
by the author m- The corresponding numerical methods were based on the 
solution of algebraic Riccati equations and certain symplectic eigenvalue prob¬ 
lems. However, computations were notoriously ill-conditioned due to singular 
Jacobians or higher multiplicity of eigenvalues. The reason for these numerical 
difficulties is an inherent rank deficiency of the positive semidefinite test matrix 
M. A rank defect one is well-known being due to the preconsistency condition, 
m.m- This rank defect increases for higher order peer methods. By an explicit 
transformation this rank deficiency is identified here exactly and an explicit 
block decomposition M = 0 0 is provided where the lower block can 
be positive definite. Positive definiteness can be checked reliably by numerical 
computations. We expect that this background may be beneficial also in the 
context of other General Linear Methods. In this paper no numerical meth¬ 
ods for finding A-stable methods are discussed, with one exception methods are 
constructed by formal computations with Maple. 

Much of the structure exploited in our analysis is due to the requirement of 
higher order. The preconsistency condition for peer methods simply reads 


Bl = l, ]1:=(1,...,1)T, 


(1.5) 


leading to one eigenvalue 1 of H. For s-stage peer methods order s — 1 is easily 
obtained for arbitrary stepsize sequences by satisfying certain matrix equations. 
For stiff problems so-called stiffly-accurate peer methods (ED containing no 
explicit terms due to the choice A = 0 have been successfully used in several 
papers, e.g. EE). For such methods and constant stepsize the conditions for 
order s — 1 correspond to the choice ([23]) 


B = {I-GE)e. 


( 1 . 6 ) 
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The matrices E and 0 correspond to polynomial differentiation and extrapola¬ 
tion at the off-step points tki = tk + hci. i = 1,..., s, k = m, With the aid 

of the Vandermonde matrix V = G the shift matrix Fq = 

and the scaling matrix D = diagj(i) they are given by if = VDFjV~^ and 
0 = VPV~^. The well-known relations P = exp(i?F(J) and 

0 = e® = exp(if) (1.7) 

between both are discrete versions of the Taylor theorem and will play a crucial 
role later on. Of course, the preconsistency condition (lESD follows from (HU) 
since V~^l = ei is the first unit vector. By choosing the nodes and the matrix G 
a stiffly-accurate peer method of order s — 1 is specified completely. On general 
grids zero stability, i.e. stability of the scheme applied to the trivial ODE y' = 0, 
may be a difficult requirement for peer methods, see |23) . For constant stepsize 
it requires that the powers of B are uniformly bounded. This leads to the 
condition g{B) < 1 where eigenvalues on the unit circle are non-defective. 

Positive definiteness of symmetric matrices is denoted hy Z y 0 or Z € 

:= {X G = V 0}, semi-definiteness by Z y 0 or Z G '.= 

{X G R®^® : X^ = X y 0}. Only real matrices will be considered but vectors 
may be complex a:: G C®, and x* = x^ denotes their conjugate transpose. In 
Section [2] we present the basic approach for establishing A-stability through a 
definiteness argument depending on an unknown weight matrix Z G 17®. It leads 
to a criterion for the generalized numerical radius of two matrices. Using a result 
of Ando an equivalent characterization by semi-definiteness of a (2s) x (2s) test 
matrix M is obtained which depends on a second matrix W G T>q. In Section [3] 
we concentrate on methods with order s — 1 and present a transformation re¬ 
vealing the rank deficiency of the test matrix explicitly. A careful analysis of the 
requirements leads to an explicit relation between the unknown weight matrices 
W and Z and an existence result for A-stable peer methods for arbitrary s G N. 
Also, the focus is moved from the existence of weight matrices to the existence 
of methods and a parametrization of A-stable methods is obtained where one 
of the weight matrices lU or Z is a parameter. The construction of practical 
methods in Sections 4 and 5 uses different techniques. For diagonally-implicit 
methods triangular decompositions and a certain triangular canonical form of 
matrices are used. The design of parallel methods with diagonal G is based on 
a representation of V~^GV from [22]. In both cases A-stable example methods 
are constructed for s = 3 and 4. 


2 A-stability and the numerical radius 

Applying the peer method (11.41) to the test equation (11.11) yields the recursion 
Ym = 'M.(z)Yra-i, where M(z) is the stability matrix of the peer method given 

by 


m(z) = (I-zG)-^(B + zA), zgC. (2.1) 


4 


Obviously, M(0) = B and stiffly accurate methods with nonsingular coefficient 
G satisfy M(oo) = 0 since ^ = 0. Hence, if such methods are A-stable they 
are also L-stable. Due to the simple product form a symmetry between the 
eigenvalues A and the parameter ^ may be used in the formulation which is well 
known from multistep methods. Obviously, for det(/ — zG) ^ 0 some A G C is 
an eigenvalue of M(z) iff 

dei{B + zA + \zG-\I)=Q. (2.2) 


This condition may be considered as an equation either for A or z with solutions 
A(z) or 2 ;(A). Hence, A-stability may be characterized in the following two 
equivalent ways where for each pair (z. A) satisfying (12.21) holds 

< 0 ^ |A| < 1, , . 

|A| > 1 ^ ?ftz> 0. ^ 


Moving A along the unit circle and solving for z gives the root-locus curves defin¬ 
ing the stability regions of multistep methods, [SI Sect.V.l]. A-stability means 
that these curves do not intersect the open left complex halfplane. Considering 
z as a solution of equation (12.21) and A as a parameter leads to the generalized 
eigenvalue problem 

z{G + ^A)x = {I - ^B)x, x^O, (2.4) 

A A 

with eigenvalue z = z(A). For A-stable methods and |A| > 1 only solutions z 
with non-negative real part are possible. This property may be implied by a 
definiteness condition with respect to a Hermitian, positive definite matrix Z. 
In the most general form of the following argument this weight matrix Z might 
depend on A, i.e. Z = Z{X) = Z{X)*. However, in this paper we consider only 
a constant, real weight matrix Z G 'D‘^ . Thus, using some nontrivial cc G C®, 
and |A| > 1, equation (12.41) is multiplied by x*[G + (1/A)A)^Z. The real part 
of the result is 

2{^z)\\Z^/^{G+\a)x\\1 =2^x*{G+^A)*Z{I-\b)x 
A A A 

=x*(G'^Z + ZG--^(A'^ZB + B'^ZA))x (2.5) 

1^1 

+ 2u\x*{ZA-G'^ZB)x. 

A 

Hence, if the right-hand side of this equation is non-negative for any A with 
|A| > 1, 5Rz must be so, too. In fact, this condition only depends on the absolute 
value = 1/|A| < 1. 

Lemma 2.1 Let Z G 'D‘^ be such that := ^{G^Z + ZG — fjL^^AJ ZB + 
BAZA)) >- Q, jjL G [0,1]. Then, a sufficient condition for A-stability of the peer 
method dn is 


\x*{ZA-G'^ZB)x\ 

x*Qfj,x 


<1 VxG e\{o}, [0,1]. 


( 2 . 6 ) 
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Proof. For a G C it holds that max{5R(j) : |A| = l//r} = f^\a\. Hence, 
requiring non-negativity of the right-hand side of (12.51) with fixed x 0 for all 
|A| = l//r means 

min (x*(G'^Z + ZG--^{A^ZB + B'^ZA))x 
/i|A|=i |Af 

+ 2^^x*{ZA-G'^ZB)x) 

A 

=x*{G'^Z + ZG- ^l'^{A^ZB + B'^ZA))x-2^l\x*{ZA-G'^ZB)x\ > 0. (2.7) 
The statement follows by considering arbitrary x ^0 and /r G [0,1] with 0: 


2/ilx*(ZA-G~^ZB)xl 
X* (G'^Z + ZG - /^^(A^ZB + BTZA)) x 


< 1 V/iG [0,1]. 


□ 

The inequality (12.61) is a condition on the numerical radius of a matrix pencil. 

Definition 2.2 Let U G and TV G (D®. The (generalized) Rayleigh quo¬ 

tient of the matrix pencil (U,N) is defined for a; G C® \ {0} by x*Ux/x*Nx. 
The numerical range (or field of values^ of the pencil is the set of all Rayleigh 
quotients, {x*Ux/x*Nx : a; G C® \ {0}}. The absolutely largest element in the 
numerical range is denoted as the numerical radius 

r([/,7V) :=max{|-^| : xGe\{0}}. (2.8) 


Remark 2.3 a) The numerieal range is a convex set containing all eigenvalues 
of the generalized eigenvalue problem Ux = XNx, 

b) For the matrix pencil U — XI the matrix N = I is omitted in the notation, 
r{U) := r(U,I^. The general form is only a simple extension of r{U) since 
any general N G 22® possesses a Cholesky decomposition N = LLA and it holds 
that r(U,N) = r{L~^UL~"^) using the shorthand notation L~^ = . The 

numerical radius r(U) has some interesting properties. It is a matrix norm 
bounded by the spectral norm, g{U) < r(U) < ||t7||2. R is not sub-multiplicative 
but it still satisfies the power inequality ^111/^112 < r{U^) < r{U)^'ik G N. 
Hence, r(U) <1 (< i) implies power boundedness (contractivity) of U, 


The dependence on the parameter p in Lemma 12.11 makes the discussion 
difficult. However, for stiffly accurate peer methods with coefficient H = 0 it 
is obvious that p, = 1 is the critical value in (12.61) . With the exception of the 
first paper on peer methods |20] only such methods have been considered in 
the literature for the solution of stiff problems due to their superior damping of 
stiff solution components. For A = 0 the condition (12.61) for A-stability can be 
re-stated with the numerical radius as 

r{G^ZB,Q) <1, Q:=^{G'^Z + ZG)>-0. (2.9) 
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If this condition is satisfied it will be with equality for any practical peer method 
due to the preconsistency condition (11.51) . In the Rayleigh quotient the precon¬ 
sistency vector 11 gives the nontrivial denominator = l^ZGl > 0 and the 

nominator ZB1\ = is the same. This critical vector 11 leads to 

tight inequalities in the new criterion for A-stability similar to those in algebraic 
stability [5]. 

The condition (12.91) is not very practical yet since the computation of the 
numerical radius is difficult. However, there exists a purely algebraic criterion 
due to Ando, see also [Ill- 

Lemma 2.4 (HI) For any matrix U G the condition r(U) < 1 is equiva¬ 
lent with the existence of a symmetric matrix W S such that 


(21-W 

I 



( 2 . 10 ) 


The statement differs slightly from the original one which used a matrix W' = 
I — W. The version (I2.10|) is more appropriate here and reveals that W itself 
is semi-definite, W G Vq. We also note that semi-definiteness in (12.101) with 
the choice W = I corresponds to the stronger condition ||f7||2 < 1- Lemma [2)4] 
yields the following purely algebraic condition for A-stability of peer methods 
with A = 0. Actually, it uses slightly weaker assumptions than Lemma ETTI 

Theorem 2.5 If there exist matrices Z and W G Vq sueh that 

+ ZG-W -G'^ZB 
-B'^ZG W 

then, the stiffly aecurate peer method disi) with A = 0 is A-stable. 

Proof. Since the matrix Q may not be definite. Lemma 12.41 cannot be used 
directly. However, considering the quadratic form of the semidefinite matrix M 
with vectors of the form , x G C®, a G C, gives 



)^ 0 , ( 2 . 11 ) 


0 <(a:*,aa;*)M 

=2x*Qx - 2ift{ax*G'^ZBx) + (|ap - l)x*Wx 

Choosing a such that |a| = 1 and 5ft(ax*G^ZHx) = \x*G~^ZBx\ the condition 
is identical with (lO) showing A-stability. □ 

Remark 2.6 For the unknown matrices Z and W the inequality (12.111) is a 
linear feasibility problem of semi-definite programming, see iHSf . For any feasible 
solution {Z, W) also all positive multiples are solutions and some normalization 
may be appropriate. 

With known matrices Z, W verification of A-stability by Theorem 12.51 requires 
only the computation of a Cholesky-type decomposition of M. Hence, for ra¬ 
tional or algebraic entries exact verification is possible. This is an advantage 
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over other criteria from the literature like the Schur criterion m requiring that 
certain rational functions be positive or bounded by one on the whole complex 
halfplane. Exact verification may be crucial since the test matrix M is rank defi¬ 
cient and semi-definiteness may not be decidable with floating point arithmetic. 


Lemma 2.7 Under the assumption cu, B1 = I, the matrix M is rank defi¬ 
cient. 

Proof. M has a nontrivial kernel, its quadratic form with (l^, vanishes. 

□ 

We note without proof that similar to m. 0 Sect.1.9] the following necessary 
conditions on W and Z may be derived, 

{B'^ - I)ZG1 = 0, Wl = B'^ZGl. 

Although the initial approach (12.71) leading to the A-stability criterion was in¬ 
dependent from the concept of algebraic stability the final criterion from The¬ 
orem [53] obtained for a constant weight matrix Z is strongly related. By using 
Sylvester’s law of inertia and special congruence transformations it was shown in 
[T5] that algebraic stability of stiffly accurate peer methods is essentially equiv¬ 
alent with the existence of a diagonal matrix D and a weight Matrix W € 
such that 

(G'^D + DG-G'^WG G'^DG-^B\ 

B'^G-^DG W 

Assuming a nonsingular coefficient G and introducing the matrix Z := G~^DG~^ 
this block matrix is congruent with M from (I2.11|) . Hence, the main difference 
between algebraic stability and the sufficient criterion for A-stability is that the 
weight matrix Z is restricted by the requirement that G^ZG = D is diagonal. 
In algebraic stability also the conditions on W and Z have changed slightly, 
W G 'D^’ and D G Vq. Hence, the results from the next sections may also be 
of interest in this larger context. In fact, by combining the congruence trans¬ 
formation from |19j with the original proof from [1] the following norm bound 
for the stability matrix M( 2 ;) (12.11) of stiffly-accurate methods is obtained. It 
verifies the assumption of the Nevanlinna theorem [Si- 

Lemma 2.8 If there exist W, Z € (D® such that M 0 in (12.111) . then for the 
stability matrix M(z) = (/ — zG)~^B it holds that 

||IT^/ 2 M(z)W^/ 2||2 < 1 Vz e c_. 


Proof. With the variables from the time step Ym — BY^-i = zGYm G C® 
applied to dm the quadratic form of the real matrix M is nonnegative. 


0<(F^,F^_i) 

=Y:,_,WYm 


(G^Z + ZG-W -G'^ZB\ { Ym\ 
'V -B'^ZG W J [Ym-J 

_i - Y*WY^ + 2^(YZG^Z(Y^ - BY^_i 
_i - Y*WY^ + 2^(z)Yf,G^ZGY^. 


)) 
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Since Sftz < 0 and ZGYm > 0 the assertion follows by = 

Y^WY^<Y^_^WY^_^ = \\W^/‘^Y^_^\\l □ 

The problem of finding suitable weight matrices D,W for (12.121) was dis¬ 
cussed by Hewitt and Hill and later by the author [TH] by considering al¬ 
gebraic Riccati equations. Hill m also used a relation with certain symplec- 
tic eigenvalue problems. However, both approaches suffer from ill-conditioning 
caused by preconsistency, see Lemma [2771 since the Riccati equation has a singu¬ 
lar Jacobian in its solutions and the symplectic eigenvalue problem has critical 
eigenvalues with higher multiplicity. Hence, numerical iterations may stall or 
suffer from severe loss of accuracy (0)- 

In this paper, we will use a different approach which explicitly reveals the 
rank defect and also changes the focus of the discussion. The algebraic criterion 
will be considered as a tool in the construction of A-stable peer methods. In 
this context, not the weight matrices W, Z of a given method are unknown in 
the first place but the whole method with coefficients B, G and the weights W, Z 
at the same time. Hence, the following discussion treats all these matrices as 
unknowns. However, additional requirements are introduced in order to narrow 
the search to methods of practical interest. 


3 Stiffly accurate peer methods of order s — 1 


For methods of higher order the coefficient matrices B and G are coupled. The 
relation for order s — 1 is given by dm), B = e- GEQ = (/- GE)e, where 0 
depends on the off-step nodes (cj) only. Thus, the A-stability condition (I2.11|) 
changes for methods of order s — 1 to 


f G'^Z + ZG-W 
-GE)'^ZG'^ 


-G'^Z{I - GE)e\ 

w ) 


>- 0 . 


(3.1) 


This condition is discussed now as a nonlinear feasibility problem for the triple 
(G, IF, Z) of unknown matrices without focusing attention on any particular 
unknown. The implicit dependence on the nodes through the Vandermonde 
matrix V is considered later. 


3.1 A rank revealing congruence transformation 

Positive semi-definiteness of the block matrix in (13.ip requires semi-definiteness 
of both diagonal blocks, e.g. IF ^ 0 in the second. Not much information can 
be gained for the first diagonal block since it has a nonlinear dependence on 
all three unknowns G, IF, Z. Now, an additional congruence transformation on 
(EH) is introduced and it is also convenient to use a transformed weight matrix 
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w := e-^we-^: 

0-^W G~^Z + ZG-W -G~^Z(I-GE)e\ / I 0\ 

1^0 e-'^J [-G'^(i-GE)~^zG~^ w y v®”^ ®”y 

{G'^ZG)E + E'^iG'^ZG) + # - W - G^Z + {G'^ZG)E 

W - ZG +E'^iG'^ZG) W 

(3.2) 

Before proceeding we note that there is the interesting class of FSAL methods 
{first same as last) [21j with a singular coefficient G. Hence, nonsingularity of 
G will be assumed only later on starting with subsection 13.31 

The transformed matrix (13.21) is better suited for the analysis than (13.11) since 
the first diagonal block depends linearly on the two matrices W and Z := G^ZG, 
but not on G alone. Furthermore, it allows to identify the rank defect of M 
accurately since the differentiation matrix E is nilpotent and also one is an s- 
fold eigenvalue of the extrapolation matrix 0. The two linear maps associated 
with the first diagonal block are considered in detail, on a matrix X G 
they act by 



Ce: X^XE + E'^X, 

Ve- X ^G'^XG- X, 0 = e^, ^ ^ ^ 

see dni. Both C and V map the set of symmetric matrices onto itself. Simul¬ 
taneous similarity transformations of the matrices E and 0 lead to congruences 
for these maps. Indeed, with X G R®^® and nonsingular (7 G M®^® we have 

U'^£e(X)E =(U~^X17)(E-^EE) + (E-^E[/)~^(U~^XU) 

=Cu-^Eu(U^XU), (3.4) 

U^Ve{X)U =Vu-^Eu{U^XU). (3.5) 

These congruence-similarity transformations will be an important tool in the 
rest of the paper establishing a similarity between different peer methods. How¬ 
ever, such methods are of practical use only if the matrices E, 0 have an inherent 
Vandermonde structure. This means that they are similar to V~^EV = E = 
DFJ and V~^GV = P by a real Vandermonde matrix V. Many proofs will be 
based on the transformed ’Nordsieck’ version with the sparse matrix E and the 
triangular Pascal matrix P = exp(i?). 

There is a strong relation between the maps £e and Ve which will be helpful. 
This relation is easily seen by considering the matrix representations of these 
maps. B^ introducing the r;ec-operator mapping matrices from R®^® to vectors 
from R® by stacking the columns, i.e. vec(X) = (xii,X 2 i, ■ ■ ■ ,Xs-i^s,Xss)~'', 
these maps may be written with the Kronecker product in the form 

vec(£E(X)) =(/ 0 E~^ + E~^ 0 I)vec(X) = Evec(X), 
vec{VE(X)) =(G^ 0 0^ — / 0 I)vec(X), 
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with the definition E := / ® + E~^ ® I. We note, that both matrices are 

nilpotent since they are similar to the strictly lower triangular matrices I 0 
E^ + E^ 0 / and 0 — I ® I. Now, since the matrices 1 0 E^ and E^ 0 1 

commute and 0 = e^, uni, the matrix of the map Ve satisfies 

0 - / 0 / 0 - / 0 / = - J 0 / 

=e^ - I = E(p{E), (3.6) 


where ip is the holomorphic function 

ip{z) = - - dt, z G C. (3.7) 

^ Jo 

The function tp is a fundamental tool in exponential time integrators, e.g. m- 
Translating (EH) back to maps on matrices we write ip(E)vec{X) = vec{^E{X)). 
Some properties of the map are proved in the following lemma directly from 
an integral representation. The lemma will use the notion of real holomorphic 
functions. This means a function q{z) with real coefficients having an absolutely 
convergent expansion in a neighborhood oi z = 0. Since such functions are 
applied to the nilpotent matrix E G only, q{E) is always a polynomial of 

degree not exceeding s — 1. 

Lemma 3.1 With the nilpotent matrix £( £ R®^® let ^e '■ R®^® —>• R®^® be the 
linear map given by 

<^e{X):= / 

Jo 

Then, with any symmetric matrix X G R®^® the following properties hold. 

a) ^E{X)E + E'^^E{X) = <d^Xid-X,i.e. 

Ve = C-E o d)£; = ^E o C-E- (3.8) 

b) $£;:P®^P®, ^E-V^o^V^o- 

c) With any nonsingular matrix U holds 

U^^e{X)U = ^u-^Eu{U^XU). 


With any real holomorphic function q satisfying q{0) = 1 holds 


q{E'^)^E{X)q{E) = <pE{q{E'^)XqiE)). 


d) For any matrix K from the kernel of Ce holds ^e{K) = K. 

e) The kernels of the maps Ce and Ve are identical. 

Proof, a) Partial integration shows that 


^e{X)E = / XEe^‘^ dt = 


JE j^^tE 


0 Jo 


^tE^Xe*^ dt 


=e^^Xe^ -X - E'^^e{X). 
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Reordering terms and recalling e-® = ©yields ^e{X)E+E^^e{X) = Ce{^e{X)) 
0^X0 — X = Ve^X). The same computations apply for ^e^XE + E^X). 

b) Definiteness follows directly from the integral representation. In fact, with 

any X € u € \ {0}, holds 

u*<i>EiX)u= [ {e*^uyX{e*^u)dt. 

Jo 

Hence, the integral is positive or nonnegative whenever X G or X G Vq, 
respectively. 

c) Let E' := U-^EU. Then, 

/*1 

U^‘^e(.X)U = / U'^ U-'^U'^ XUU-^y^Udt 

y<.E'y u'^ XU dt. 

A matrix of the form U = q{E) commutes with E leading to E' = E. 

d) Let K be any matrix satisfying KE = —E^K. Then KE^ = {—E~^yK, j G 
N, and from the series expansion of the exponential function follows ATe*^ = 

K_ Hence, 

^e{K)= j U^^KU^dt=f = K. 

Jo Jo 

e) The matrix representation of the linear map $_e is (/^(E). By assertion c) we 

may consider the map with the transformed matrix E~^ 0 / + / ® E^ which is 
strictly lower triangular. Hence, ip{E'^ ® ® E~^) is lower triangular with 

unit diagonal and, hence, nonsingular. By (13.81) ^e^Ue and Ce have the same 
kernel. □ 

Remark 3.2 Since E is nilpotent the map d>£; = ^^sily eval¬ 
uated in practice. For X G the recursion Xq := X, 

Xk := CE{Xk-i) = Xk-iE + fe = 1, 2,... 

terminates with X 2 S -1 = 0. This follows inductively from the fact that each 
application of the transformed map Cp introduces one more trivial antidiagonal 
(southwest to northeast) in Xk starting with = 0. Hence, 

2s-2 ^ 

= E 7i-rTiT^‘- <*■“> 

k=0 ' '' 

With the definition (13.31) at hand the leading diagonal block Mu of the test 
matrix dUl) may be written in a more compact form. Its semi-definiteness is 
necessary for the semi-definiteness of the whole test matrix itself. With the 
factorization (j3.8l) this requires that 

Mil = Ce{G^ZG) -Ve{W) = Ce{G^ZG - <^e{W)) E 0. (3.10) 
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It is important to note that due to the singularity of the map Ce semi-definiteness 
of an image /is(X) ^ 0 is a strong restriction for both X and the image Ce{X). 
These conditions follow easily from the simple structure of the transformed ma¬ 
trix E = DFJ. Here, the image U = £^{X) is given by 

Uij = {j - + {i - l)xi-ij, l<i,j<s, (3-11) 

where elements with an index zero are missing. Hence, the first diagonal element 
of the image is always zero, un = 0 , and the first row and column must vanish 
completely for U to be semi-definite. In fact, only a block on the southeast corner 
of U may be nontrivial. The precise structure is given in the next lemma. 

Lemma 3.3 Let X G be symmetric and let U := Cj^{X) > t). Then, 

Uij = 0 for min{i, j} < [(s -I- I)/2J and also Xij = 0 for i + j < s. 

Proof. By induction over the rows i = 1,... ,k := [(s -|- 1)/2J we show that 
= 0; J = b ■ • ■, s and Xij = 0, j = i,... ,s — i. If some diagonal element 
oi U G T>q vanishes, uu = 0, then by semi-definiteness all elements of row and 
column number i vanish, as well. Hence, the first row and column of U are zero, 
since un = 0. In (13.111) this gives 0 = uij = (j — l)a;ij_i, j = 2,..., s, i.e. the 
assertion for i = \. Now, by induction with i>2 we get uu = 2(i — \)xi-i^i = 0 
for 2i — 1 < s, i.e. i < k, and this leads again to Uij = 0, j = i,..., s. In 
(|3.1ip follows again 0 = Uij = (i — l)xi-ij + (j — l)xij-i = (j — l)xij-i for 
j=i,...,s + l — i leading to xu = 0 for i < i < s — i. □ 

It is also possible to describe the kernel of the map on the set of symmetric 
matrices explicitly. By (13.1111 and also since ejE = 0, the rank-one matrix esej 
is always an element of this kernel. Other elements can be derived from the 
proof of the last lemma. As in the base case of induction the first row of 
any matrix K satisfying Ce{.K) = 0 vanishes with the possible exception of 
the last element, fey = 0, j = 1,..., s — 1. Now, (13.111) is a linear one-step 
recursion along antidiagonals and it follows that also kij = 0 for i + j < s. 
On antidiagonals with i + j > s, kis 7 ^ 0, this recursion leads to elements with 
alternating signs. Hence, corresponding solutions K are symmetric only if the 
antidiagonal contains a diagonal element, i.e. ii i + j is even. This result is 
formulated as a simple lemma which by dSH) remains valid with any matrix E 
similar to E. 

Lemma 3.4 The kernel of the map Ce on the set of symmetric matrices has 
dimension [(s-|-l)/ 2 j. 

A discussion is similar for the pre-image of semi-definite matrices U. Since 
the map Vp, is sometimes more convenient in the construction of methods (see 
Section [5]) the following result is formulated for this map. 

Example 3.5 For s = 4 the general symmetric pre-image of a semi-definite 
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matrix U = Vj^{X) ^ Vq is given by 


X = 


/O 0 
0 0 
0 + 


0 

i'«33 

2^33 

* 


-|U33 \ 

— fa;33 — 1^33 + 1^34 

3^33 — |m34 + 1^44 

2^44 / 


(3.12) 


This means that only it 33 ,U 34 ,M 44 of the image may be nontrivial and X 33 and 
0:44 are parameters of the kernel oiVj^. 


3.2 Explicit feasible solutions 

Returning to the feasibility problem the test matrix is rewritten with 

Lemma ED in the form 


M = 


(Mu 

\MJ2 


Mi2\ 

W j 


( Ce{G'^ZG-<^e{W)) W - G'^Z + {G'^ZG)E 
\W - ZG + E'^{G'^ZG) W 


(3.13) 


By Lemma [3731 Ml 1 has a rank defect [(s + 1)/2J and by semi-definiteness M 12 
must have the same, at least. Equation (I3.13|) has been used as a definition of 
the blocks My so far. However, with a matrix M ^ 0 having an appropriate 
kernel it may be considered as a set of equations where M is the slack variable 
of the feasibility problem. Peer methods are suited for stiff problems only if 
the eigenvalues of the coefficient matrix G have nonnegative real part. This 
property follows directly from (13.131) . 

Lemma 3.6 Let G G and Z G , W GVq (or W G ) satisfy (I3.13|) 

with fS/I G T’o*- Then, all eigenvalues of G have nonnegative (positive) real part. 

Proof. The equation for the off-diagonal block is written as G^Z = {G~^ZG)E+ 
W — Mi 2 . Adding its transpose leads to 


ZG + G~^Z =(G^ZG)E + E'^iG^ZG) + 2 W - M12 - 
=©■^1X0 -k Mil - M 12 - A/T -k W 
^©■^IT© ^ 0. 


The last step used the fact that for M G also Mu — M 12 — MJ 2 EW GV^. A 
congruence transformation with the inverse factor of a Cholesky decomposition 
Z = hzL^z IfiS-ds to the symmetric matrix 

lIGLz'^ -k L-/g'^Lz > > 0 . 

Hence, the symmetric part of the matrix is positive semi-definite and 

all its eigenvalues have nonnegative real part. Of course, is similar to 

G. The argument for W GV^ is the same. □ 
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The leading blocks Mu and M 12 are highly rank deficient by Lemma [3.31 
Hence, simple solutions may be obtained where these two blocks are chosen zero 
to start with. This means that solutions Z £ W £ G G of the 

following system are sought: 

Le{G'^ZG-^e{W))= 0, 

W - G'^Z + {G'^ZG)E = 0. ^ ’ 

The second equation is not solved for W to eliminate it from the system since 
symmetry and (semi-) definiteness of W are crucial. Enforcing these properties 
would lead to difficult side conditions. 

The first equation of the system (13.141) contains the unknown matrix Z = 
G^ ZG. For nonsingular G this matrix is definite again and may be considered to 
be the unknown weight matrix itself. Then, the second equation takes the form 
0 = W — Z(G~^ — E). The first equation in (13.141) means that Z — $£;(IT) is an 
element of the kernel oi Ce, for instance zero. Now, choosing W 0 and solving 
the first equation with Z = ^e{W) and the second for H := G~^ = E + Z~^W 
yields a special feasible solution. This gives the following general theorem. 

Theorem 3.7 For each s £ N there exist stijfly accurate A-stable peer methods 
with s stages and order s — 1. 

Proof. Let s £ N and let c,, i = 1,..., s, be an arbitrary set of mutually 
distinct nodes and V its Vandermonde matrix. Then, with E := VDEqV~^ 
and any matrix W € also Z := ^e{W) £ 2?® by Lemma [3.11 and the pair 
{W,Z) solves the first equation in (13.141) . Since the inverse of Z exists by 
definiteness, a solution of the second equation \s El = E + Z~^W. With these 
matrices the test matrix (13.21) is semi-definite with one nontrivial block IT £ 22® 
and the matrix 22 is nonsingular by Lemma 13.61 Hence, the matrix G = 
and B = {I — GE)Q from (11.61) define an A-stable peer method of order s — 1 
according to Theorem 12.51 and (13.21) . □ 

Unfortunately, Theorem 13.71 has limited practical relevance since the coeffi¬ 
cient matrix G will be a full matrix, in general, leading to prohibitively expensive 
peer methods. Turning attention to diagonally implicit or parallel methods the 
discussion is restricted to the simpler case of a nonsingular coefficient matrix G. 


3.3 Methods with nonsingular coefficient matrix G 

For nonsingular G it is convenient to rewrite the test matrix ()3.13l) in terms 
oi Z = G^ ZG £ 22® and El = G~^. Now, the test matrix has a generic form 
depending on the quadruple (2?, 22, W, Z). This form remains unchanged under 
congruence-similarity transformations and is denoted by 


MiE, H,W,Z) 


/£e{Z)-Ve{W) W-Z{H-E) 
\W-{H-EyZ W 


(3.15) 
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Remark 3.8 Congruence-similarity transformations with block diagonal matri¬ 
ces diag(U,U) establish an equivalence relation between test matrices of the form 
(|3.15l) . In fact by (13.41) , (13.51) it follows that 

M{E,H,W,Z)hO ^ M{U-^EU,U-^HU,U^WU,U^ZU)>Q (3.16) 

for nonsingular U £ This property will be helpful in the construction 

of efficient peer methods where H = G~^ has lower triangular or diagonal 
form. However, the resulting methods only have high order if their E-matrix 
is Vandermonde-similar to E = DEJ . This distinguished Nordsieck version 
using U = V is used frequently and the actual value of A4 is denoted as 
M = A4{DFJ , H ,W, Z). In fact, the rank defect of the first block has been 
identified by Lemma \3.S\ for Mu. By semi-definiteness the same rows as in Mu 
must also vanish in M 12 leading to an explicit decomposition in block diagonal 
form M = OfeSM^i with k = [(s + l)/2j mentioned in the introduction. It will be 
seen in some examples that the nontrivial block may be definite, Md £ 2 ?^®“^, 
which can be checked reliably with floating point arithmetic. The construction 
of A-stable methods will make use of this explicit block decomposition. 

With (j3.15|) the original feasibility problem is again written in the form 

= (3,17) 

with the matrix M as slack variable. With the results of the last subsections it is 
possible to derive a rather explicit parametrization of solutions of the feasibility 
problem Ai(E,H,W, Z) >: 0. We note that in the Nordsieck form pre-images 
X with -Pe {X) = Mil £ Rq computed easily, see Example 13.51 

Lemma 3.9 Let E £ R®^® be nilpotent and Wq £ R®. With a matrix K = 
from the kernel, Ce{K) = 0, and N = N~^ £ K®^® such that Ve{N) = Mu £ 
££;(Rg) (~1 Vq let W := Wo — K — N € Then, for any matrix M 12 £ K®^® 
such that M £ Rq® the triple of matrices 

Z = ^e{Wo), W = Wo- K - N, H = E + Z-\W + M 12 ), (3.18) 

is a feasible solution M{E, H, W, Z) ^ 0 with Z y 0. 

Proof. By Lemma 13.11 the equation for the first block may be replaced by 
Ce{Z) — Ve{W) — Mil = P,e{Z — ^e{W -\- Nf) = 0 and its general solution 
is Z = ^EiW -\- N) -\- K = $e(W -\- N -\- K) with Ce{K) = 0. Assuming 
Wo '.= W -\- N -\- K & T>^ gives Z £ R® by part b) of Lemma 13.11 Now, since Z 
is nonsingular, the off-diagonal block equation may be solved uniquely for H. 
□ 

Before turning to the construction of practical methods with restrictions on the 
shape of H = G~^ two interesting cases of transformations (13.161) of known 
solutions are considered. First, part c) of Lemma 13.II shows that transformation 
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matrices commuting with E lead to simple relations. In fact, congruences with 
a real function q{E) leave the kernel invariant: if Ce{K) = 0, then 


CE{q{E'^)Kq{E)) =q{E'^)Kq{E)E + E'^ q{E'^)Kq{E) 
=q{E'^){KE + E'^K)q{E) = 0. 


Hence, any solution (13.181) without slack elements (Mii = 0, M 12 = 0) yields a 
new solution by the similarity transformatioir 

q{E)-^Hq{E) = E+{<^>E{q{E'^)Woq{E))y\q{E'^)Woq{E)+K). (3.19) 

Hence, if E is already Vandermonde-similar to .E by ill = V, (I3.19|) 
may be used as a post-processing to change the shape of ii. Effects of this 
transformation on slack terms have to be discussed on a case-by-case basis. 

Of fundamental interest may be an exceptional transformation eliminating 
the weight matrix Z. It yields a simple, distinguished representative of the 
whole equivalence class of relation (13.161) . It uses the Cholesky decomposition 
Z = LzL^z introduces the transformed matrices W := L^W, E := 
L^ELy, H ■= L\E[L~y. Then, the matrix M from (j3.17|) is congruent with 


MiE,H,I, W) 


(Cyi)-VE{W) W-H + E\^_^ 

\W -H^ + E'^ W J ■ 


(3.20) 


depending linearly on W and ff. Since Z has been eliminated an explicit so¬ 
lution representation for W is required now. This is obtained with the inverse 
map of ^e satisfying $£;(dt£;(V)) = X for X e T>^ and nilpotent E. Similar 
to (13.61) its matrix representation is given by 'ip(E)vec(X) = vec(4'£(X)) where 
iIj is the reciprocal of the function ip from (ITTI) . 


V’(z) 


1 

ip{z) 


z 

2 


z /e^ + l\ 

2 V - 1 / 


2 

2 


+ i-^(-i)'= 


h 2k 

(2fc)! 


(3.21) 


The series expansion of ^ has only one odd term —z/2 and the even terms 
define the Bernoulli numbers /?fe, and the radius of convergence is tt. In fact 
the expansion corresponds to the Magnus series ip = dexp“^, see [SJ §IV.7]. 
Using e(^) representing the general feasible solution in (13.201) without 
slack, i.e. Mu = M 12 = 0, shows that the skew-symmetric part of the matrix 
H is identical with that of E. 


Lemma 3.10 Let E G be nilpotent and W = — K G 21® with 

K = and Cg{K) = 0. Then, a solution with Mu = M 12 = 0 in (I3.20p is 
given by W and H = E + W = Hsym + Hskew, where 

Hsym =yE + Ey + + K = I- + k (3.22) 

Hskew —Hskew — H )• 
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Proof. By Lemma [XT] holds C^{I) = i;{I)) and the 

general solution of £^{1) — Vj^iW) = Pb('I'^(/) — W) = 0 is VL = + K 

with a symmetric matrix K from the kernel £^{K) = Vj^{K) = 0. Solving 
Mi 2 = 0byiL = i^ + 'I'^(/) the linear term from the expansion (13.211) combines 
with the term E and gives the skew symmetric matrix E— ^£j^{I) = ^{E — E~^). 
Since the remaining terms are symmetric the assertion follows. □ 

The lemma shows that introducing the map 'I 'e provides additional insight in 
the structure of solutions. However, since T e does not preserve definiteness the 
assumption W had to be added explicitly. Using the definiteness preserv¬ 
ing map $E was more convenient for the theoretical result of Theorem 13.71 
After choosing the matrix E Lemma [3.101 provides a simple feasible solution 
which may be complemented by slack terms Mij ^ 0 , j = 1 , 2 , having the 
proper rank. From it the whole equivalence class of A-stable methods can be 
obtained by similarity (13.1611 . However, it is still necessary to identify methods 
of higher order having a Vandermonde structure with E = VDEJV~^ and being 
of practical interest with a triangular or diagonal matrix H = G~^. For these 
two cases different approaches are used. 


4 Diagonally implicit peer methods 


The methods from Theorem are not suitable for practical computations 
since G G is a full matrix, in general. In this case time steps of (HJl are 

too expensive requiring the solution of a nonlinear system of s • n simultaneous 
equations. For methods with a lower triangular coefficient matrix G the s stage 
systems are decoupled. Such peer methods have been discussed in [2], [2T] . 
Singly-implicit methods with constant main diagonal = 7 , i = 1,..., s, have 
the additional benefit that costly matrix decompositions may be computed only 
once per time step. 

Similarity transformations with lower triangular matrices do not destroy the 
lower triangular form of H. A first transformation uses the L-factor of an explicit 
LU decomposition V = LyUv of the Vandermonde matrix V [24] given by 


C/7 


{Ly^) 


'ij 


/I -Cl C1C2 -C1C1C2 ...\ 

1 -(C1-I-C2) C1C2 -I- C1C3 -I- C2C3 ... 

1 -(ci -I-C 2 -I-C 3 ) 

V / 


n 

fc=i 





(4.1) 


(4.2) 


In fact, this is the matrix version of Newton interpolation with the nodes (c^) 
since L~^ contains the factors for the divided differences and U~^ contains the 
elementary symmetric functions. The factor Ly does not have unit diagonals, 
all its entries depend on the differences between nodes only. 
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It is possible to obtain a final solution H with lower triangular form obeying 
(|3.18|) with E = VEV~^ from a simpler, equivalent version. In fact, a similarity 
transformation with the factor Ly yields 

Ly^HLv =Ly^ELv + {<^L-^ELviLlWoLv))~"Ll{W + Mi2)Lv, (4.3) 

where Ly^HLy also has lower triangular form while Ly^ELy = UyDFjUy^ = 
UyEUy^ is strictly upper triangular. In fact, its entries depend only on differ¬ 
ences of the nodes and for s < 4 this matrix is explicitly given by the leading 
submatrices of 


Ly^ELy = UyDFjUy^ 


/o 

V 


1 Cl - C2 
0 2 
0 


(ci - C 2 )(ci - C3)\ 

Cl + C2 - 2 c3 

3 

0 / 


(4.4) 


We note that the terms JtT-l-TV in the matrix W, see ()3.18l) . contribute to (14.31) a 
matrix Lj^{K + N)Ly = Ry'^{K + N)Ry^ where X = K N is given explicitly 
by (13.121) for s = 4. Since upper triangular congruence transformations do not 
change the zero pattern in (13.121) the kernel and slack terms Ly{K + N)Ly in 
(j4.3p only contribute to the southeast antidiagonals. 

It is clear that one can go back from a lower triangular solution Ly^ElLy of 
the form (14.31) with a suitable strictly upper triangular matrix (14.41) to a solution 
()3.18l) by the similarity transformation with the factor Ly of the Vandermonde 
matrix. The open question remaining is to find solutions such that the if-matrix 
has the form from (14.41) . UyDFjUy^ = UyEUy^. Here, the following lemma 
will be applied to the Z-free versions (I3.20p . (I3.22p . It is proved only for the 
non-defective case. 


Lemma 4.1 Every non-defective matrix A £ possesses a triangular canon¬ 
ical form 


A = UaLaU^^ (4.5) 

where La is lower triangular and Ua is upper triangular and nonsingular. 

Proof. Let A = XJX~^ be the Jordan canonical form of A where J is diagonal 
by assumption. Then, with some permutation matrix H there exists the LU- 
decomposition = IlLxUx and the factorization A = JIlLx)Ux 

proves the assertion with La = Lffn^JIlLx and Ua = Ux\ □ 

The canonical form may not exist in some exceptional cases like a Jordan 

block , 6 yf 0, in dimension 2. In fact, the lemma is related to the canon¬ 

ical form computed by the classical LR-iteration of Rutishauser [26] . However, 
the existence of this decomposition can be used as a criterion in the search for 
methods. Obviously, the main diagonal of La contains the eigenvalues of A. 
The diagonal elements of Ua may be chosen arbitrarily but non-zero adding 
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some degrees of freedom in the forthcoming application. This triangular canon¬ 
ical form will be used to construct or reconstruct A-stable peer methods from a 
Z-free representation (see Lemma rd.lOl) and the knowledge of the matrix E and 
nontrivial kernel and slack terms alone. We consider this approach primarily as 
a theoretical tool using formula manipulation or exact arithmetic since numer¬ 
ical computations may be too inaccurate for multiple eigenvalues. The input 
of the following algorithm is the strictly upper triangular matrix E and kernel 
and slack elements from M\i which are combined in a matrix X. Of course, 
the choice of such terms is restricted by the definiteness oi W = — X. 

The image V^{X) is required to be semidefinite and has the zero structure of 
Lemma 13.31 The two nodes ci, Cg may be chosen freely, where the choice = 1 
is convenient in practice. 

Algorithm. 

Input: Strictly upper triangular matrix E € ei,i+i ^ 0 , i = l,...,s- 1 

and X = X~^ such that Vp;{X) G Vq. 

Steps: 

1. Compute W = '$— X and H = E + W, check definiteness VL £ 2?®. 

2. Compute the triangular canonical form (14.51) for H = UhLhU]j^ with 
the diagonals of Uh chosen as uu = {i — i)!/n ;=2 Cj-ij, starting with 
Mil = 1- 

3. Apply the congruence-similarity transformation with Uh to (j3.22p by 
computing E' = UJj^EUh, W = UJjWUh, Z' = UJjUh satisfying 
Lh = E' + {Z'')~^W'. The matrix H' = Lh is known from step (2). 

4. Identify E' with (j4.4l) . i.e. solve the equation E' = UvDFjUy^ for the 
differences of the nodes Ci — ci, i = 2,..., s — 1. In the first superdiagonal 
the identity eC_|_i = i, 1 < i < s is satisfied by the choice of the diagonals 
of Uh in step ( 2 ). 

5. Apply a final congruence-similarity transformation with the Vandermonde 
L-factor Ly from (14.21) to the matrices in order to obtain the coefficients 
of the A-stable method 

H = LyH'Ly\ E = LyE'Ly^ = VEV-^, 

and the weight matrices W = Ly'^W'Ly^, Z = Ly'^Z'Ly^. 

The procedure is based on the following congruences of the equation H = E+W : 

H = UhLhU]j^ = E -h W 

{UIUh) Lh = [UIUh)Uh^EUh + {UIWUh) 

II 

Z' H' = Z' UvEUy^ + W' 

{Ly^Z'Ly^){LvH'Ly^) = {Ly^Z'Ly^){LvUvEUy^Ly^) + {Ly^W'Ly^) 
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Remark 4.2 There are two critical points in this algorithm. For methods of 
practical interest eigenvalues have to he real. Complex eigenvalues can easily 
he avoided in Step 2 with the criterion of Sylvester IZ3 Th. 1 . 4 . 4 ] by computing 
traces of powers of H and a LU decomposition of a certain Hankel matrix. More 
crucial is Step 4 since (gH) has only s — 2 degrees of freedom for equating the 
(* 2 ^) elements e'^-, i + 2 < j < s. Hence, for s > 4 the algorithm will only 
he successful if the input E satisfies certain unknown side conditions. However, 
the algorithm may he used to decode the full triple {H, W, Z) plus nodes of a 
known A-stable method from a compact representation consisting of the trian¬ 
gular matrix E and optional kernel or slack contributions. 


Only for three-stage methods Step 4 consists of one single equation which can 
always be solved by C 2 = ci — This is demonstrated now by constructing 
an A-stable singly-implicit three-stage method. 

Example 4.3 For s = 3 the algorithm is used with E and one element ke^ej, 
it > 0, from the kernel of This ansatz gives R = ^{E — E~^)1— 
-p ke^ej. The conditions for a threefold eigenvalue rj oi H can be 
solved for the parameters k, 613 and 623 with expressions containing square roots. 
Enforcing positive radicands and positive diagonals of the LU-decomposition 
of W leaves nontrivial admissible regions in the ( 77 , ei 2 )-plane. Looking for 
parameters with small radicands leads to the choice rj = 5/2, ei 2 =3/2 yielding 
k = 174/175, ei 3 = -^v^, 623 = and 



The factors of the triangular canonical form (14.51) of H are 


Lh 


( ^ \ 

-81-9\/65 5 

40 2 

27(4+^65) 9(\/65-17) 5 / 

\ 70 28 2/ 


Uh 


V 


11-V65 

3 


ll-\/65 \ 

f 
18 
V35 
18 / 


The diagonal elements of Uh are chosen to give = i,i = 1,2, for E' = 

Ufj^EUn and equating the only general element e]^ = ci — C 2 from (14.41) gives 
C 2 = Cl -I- (259 — 21-\/65)/84. Fixing the free nodes as ci = —1, C 3 = 1 the 
final transformation yields H = LyH'Ly^ and the methods parameter matrix 
is obtained as its inverse 


G = 


( 5 

207+15^65 

500 

I -3672+1800^65 
\ 30625 


0 


2 

5 

-468+180v^ 

1225 
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The transformed weight matrix Z = LyZ'Ly = {ULyYULy is given by 


Z = 


/ 129447+7757V65 

345744 

* 


124345-9465V64 

749112 

760215-18615 V65 

1997632 


V 


* 


-3315-41v^ \ 
91728 1 

-3665+561V65 

122304_ 


339 + 29 \^ 

7488 


and W = LyW'Lv = {ULyYWULy can be computed by IT = ^e{Z) + 
^here k = (841 + 87V^)/49920 and = ((13 + 9v%5)/98,-(111 + 
9-\/^)/98,1) is a multiple of the last row of Ly. With these data semi-definiteness 
of the test matrix M (I3.17p can be checked exactly. Actually, the test matrix is 
zero with the exception of the last block M 22 = IT by construction. 

It is difficult to apply this construction with s > 4. First, the occurring 
algebraic equations can no longer be solved with simple square roots. Secondly, 
the matrix E' contains 3 general elements while (14.41) has only 2 degrees of 
freedom. In order to show that singly-implicit 4-stage A-stable peer methods 
exist we give an example computed by some numerical iteration in the Nordsieck 
form. It uses a slack matrix Mu = diag(0,0,1,2) of rank two and kernel 
46464 . In this form numerical verification of A-stability is reliable since the 
nontrivial part of the test matrix M is definite. In order to save space the 
Z-free formulation of this method is given. 

Example 4.4 An A-stable singly-implicit 4-stage peer method may be recon¬ 
structed with the algorithm from its compact form W ='^f^{I) + X with 


E = 


(0 

V 


1.828746674 -2.100327482 
0 1.159613679 

0 


-0.0374802335\ 

0.0397218240 

1.105306335 

0 


In this form the nontrivial part of the slack matrix T’e{X) is 


/ 1.124278133 -0.2749000029\ 

1^-0.2749000029 0.3724460117 ) 

and the kernel is determined by X 33 = —0.845690710, 0:44 = —0.7361686650. As 
mentioned above the fourfold eigenvalue rj oi E[ and H can not be computed with 
sufficient accuracy, its accurate value is 77 = 1.80350113085004, and the nodes 
are (c,) = (-0.889874593986289,0.522100340305431,-0.297184898847891,1). 
We do not go into more detail here since this method may not be of practical 
interest because efficient methods should satisfy several criteria like those listed 
in [21]. The construction of efficient A-stable peer methods will be the topic of 
a forthcoming paper. 

In principle, the feasibility problem for diagonally implicit methods contains 
so many parameters that one may expect that the algebraic criterion for A- 
stability can be satisfied for larger stage numbers s as well. 
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5 Parallel peer methods 


Due to an explicit representation of the matrix G the design of parallel meth¬ 
ods can be completely discussed in the Nordsieck version M{DFq , H ,W, Z). 
Although A-stability is discussed only for constant stepsizes here, peer methods 
should perform well on general time grids. Then, zero stability of stiffly-accurate 
parallel methods requires that the diagonal matrix G has non-constant diagonal 
elements and does not commute with V, see [53]. Using the Nordsieck version 
it is necessary to identify matrices G = V~^GV which are similar to a diagonal 
matrix G by a Vandermonde transformation V. The following necessary con¬ 
dition is related to a representation from [^. It uses the standard notation of 
square brackets for the commutator of matrices. 

Lemma 5.1 Let V £ be a Vandermonde matrix with distinct nodes ci,..., 
and let p{x) = {x — ci) ■ ■ ■ {x — Cg) = a:® -I- node polynomial. 

The vector p = (pi) £ M® contains its coefficients. //G £ R®^® is diagonal, then 


G = V ^GV satisfies 


[Uo-peJ,G] =0. 


(5.1) 


From (EH) follows the rank condition 



(5.2) 


on the submatrix of dimension (s -|- 1 ) x (s — 1 ) indicated by the subscripts. 

Proof. The diagonal elements of G may be considered as function values of 
some polynomial g of degree s — 1 at the nodes, i.e. gu = g(ci), i = 1 ,..., s. 
Introducing the matrix G = diag(ci) this means that G = g{C). Since the nodes 
are the zeros of p it is well known that CV = VF holds with the Frobenius 
companion matrix F ■.= Fq — pej. Hence, V~^GV = V~^g{G)V = g{F) is a 
function of F and commutes with it. Writing (15.ip as [Fq, G] = pejG — GpeJ 
shows that the rows of [TbiG] are multiples of the last row G itself with the 
exception of the last column. □ 

With given node polynomial the condition (15.11) is a simple linear system. 
The second version (15.21) . however, allows to characterize candidates for trans¬ 
formed diagonal coefficient matrices G without knowledge of the node polyno¬ 
mial. Rank 1 can be enforced by a set of 2 x 2-determinants leading to quadratic 
restrictions on G. The coefficients of the polynomial can be recovered from the 
equations ([Fo,G] — pejG)i...s-i = 0 . Actually, if Cs = 1 is one of the nodes 
which is convenient in practice, then l^p -f-1 = 0 and the matrix from (15.2p has 
trivial column sums, (]l^[Fo, G] -I- ejG)i...s_i = 0. 

Example 5.2 Construction of a parallel 3-stage peer method. The essential 
parameter of the construction is the matrix W £ and the ansatz 



Kl L) 

2ki 0 ^3 


2h;i 

0 
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Dw = diag(l, \d 22 , and q{z) = I + qiz + q 2 z‘^ is used, i.e. Mu = 0. It 

takes advantage of property (13.191) and leads to simpler equations in the side 
conditions than a general Wq € (D®. Condition (15.2|) is enforced for the inverse 
H = G~^ by the 3 determinants 


:= 


[Fo,HU 

hsi 


[Fo,Hh 

^32 


i = 1,2,3. 


The choice ki = 1/3 simplifies the conditions and ^22 := Qqf — 12(72 — 1 cancels 

both (5i and 62 - Solving (^3 = 0 for K 3 and annihilating the sum in the first 
column (]l^[^o, ff] + = 0 with c ?33 leaves some region in the (( 71 ,( 72 )- 

plane where Z = $g(Wo) and W are definite. Finally, looking for rational 
solutions with small denominators suggests the choice 91 = 3 and 52 = 17/4. 
The resulting nodes are (cj) = (0,2,1) and transforming back the coefficient 


form (13.171) are 


^=12 


G = 

diag(|: 

20 5 '1 

’ 29’ 111 

23 

20 

-50\ 

20 

37 

-52 , 

-50 

-52 

116 / 


Z = 


5 _5N 

3 2 

5 -5 

-"5 _ c 20 
' 2^3 


So far, the system (I3.17|) has been discussed as a nonlinear problem for H. 
However, if a peer method of order s — 1 is given by fixing the nodes and the 
matrix G or H = G~^ , (13.171) is a linear system for the unknown weight matrices 
Z and W where only positive definite solutions are of interest. This corresponds 
to the original approach to the feasibility problem from the literature, e.g. [9]. 

A simple class of implicit parallel methods has been discussed in [23] with 
the focus on zero stability for non-uniform grids. With the diagonal matrix G 
of the special form G = gol + giC, go,gi € M+, G = diag(ci), the transformed 
matrix B is upper triangular leading to a simple proof for uniform zero stability 
if the parameter 771 lies in a certain interval containing the choice 771 = 2/5 
for s = 4. We now verify A-stability of a 4-stage method of this form with 
prescribed nodes. 

Example 5.3 Construction of a parallel 4-stage method from the class dis¬ 
cussed in [53]. In order to obtain rational solutions the parameters from 
are approximated by 771 = 2/5 and nodes c = (—1,—2/5, 2/5,1). Only 770 is 
considered a variable in the design. With Z = {zij) € Zn = 1, the general 
ansatz W = ^ j^{Z) + X is used with X = K + N from (13.121) . Lemma (SHI shows 
that for s = 4 the first two rows and columns of Mu vanish. This restriction 
is obeyed already by using the matrix X in the form (13.121) . However, semidefi¬ 
niteness also requires that the first two rows of M 12 be zero. This gives a linear 
system with 8 equations which can be solved for the parameters 5 i 2 , ^ 13 ,..., Z 24 
and M 33 ,U 34 from (I3.12|) . Due to the upper triangular form of E the number of 
parameters in the entries of M = (rhij) increases from top-left to bottom-right. 
For this reason the parameter U 33 is required since it is the only free parame¬ 
ter in the equation rh^s = 0 besides go. These dependencies also facilitate the 
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choice of the remaining parameters by enforcing non-negativity of the minors of 
M step-wise while keeping a nontrivial admissible interval for the parameter go. 
At the end it was found that the LU-decomposition of M indeed has 6 positive 
diagonal elements for go G [0.789,0.822]. Hence, with go = 4/5 the method is 
proven A-stable and its coefficient matrix is G = diag(2/5,16/25, 24/25, 6/5). 
Since the weight matrices have smaller numbers in the Nordsieck form we show 
this version. 


/I 

0 

0 

5 

6 

7 

6 

1429 

2250 

473 \ 

2250 

,W = 

* 

1 

2 

1 

-1 

19 

12 

^ \ 
12 \ 
679 

4500 

* 

* 

5 

6 

* 

* 

8329 

1125 

833 

100 

\* 

* 

* 

16 yi 


\* 

* 

* 

69 , 
2 / 


This example may provide information of broader interest. It shows that it 
may be possible to construct feasible weight matrices for parallel methods with 
four and maybe more stages even if the method is specified almost completely. 
This indicates that the sufficient algebraic condition for A-stability from The¬ 
orem [23] may be rather sharp already with a constant, A-independent weight 
matrix Z as used in (12.51) . Also, the test matrix M has rank 6 which is maximal 
according to Lemma 13.41 for s = 4. Rank defect two means that the Nord¬ 
sieck version has an explicit representation M = Ofc 0 with k = [(s + 1)/2J 
where the nontrivial part is positive definite. Definiteness of may be verified 
reliably in floating-point arithmetic. 


A Detailed computations 


Space was saved in the 4 examples of the preceding sections by presenting the 
essential data of the constructed methods only. In order to facilitate verihcation 
of the criterion the different transformations for all examples are reproduced in 
detail now. 

Details for Example 14.31 The triangular coefficient matrix G and the weight 
matrix Z have already been shown, the second weight matrix for (13.151) is 



/ 76249 1 

63211 

33793 

33073 

541 1 2945 , 

w = 

/ 76832 ' 

691488^ 

115248 
1887495 1 

90585 

4704 ' 183456^ 
16225 661 

1 

* 

3995264 ' 

3995264 ^ 

* 

244608 81536 

121 , 37 3 

4992 ~ 14976 ^ 

and the 

coefficient matrix B is 

given by 




/ 817 

ggOoV^ 

+ 1 1 

675 

A _ 

80 208 '' 

_ J89. , 

2000 ^ 5200 '' 

649 1 63 

B 

/ 1960 

_ I 1371 

1 6125 

\ 55737 

10192^ 

1917 

50960 '' 

12393 


\ 245000 

19600 

254800^ 

10000 ' 1040^ 


The eigenvalues of B are = I, 0.4569, 0.0455 in accordance with preconsistency 
and zero stability. By construction the test matrix M in (I3.I5D has three zero 
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blocks and only M 22 = kk 0 is definite. Going back to the original formulation 
from Theorem 12.51 the transformed weight matrices are H^ZH = 


Z = 


/ 3389263 
768320 


I 


and ©■^ire = 


3518461 

6914880 

* 


- 


W = 


1163 

10976 


767 

98784 

* 




902677 

577413 

460992 

1997632 

18279351 

1 1358121 

7990528 

' 7990528 


* 

7585 

1647 

32928 


422145 1 

4383 

570752 ' 

570752^ 


* 






265 I 1415 
672 26208 

34595 _ 1215 


2825 I 
9984 29952 


23296 

^^65 


21 

7219 


\ 

1_1648 'V^ 1 • 


34944 _ 

1009 I 253 . / 

4992 ' 14Q7R VniJ / 


14976 


The smallest eigenvalues of Z and W are = 0.204 and = 0.0167. The original 
test matrix M (12.111) is a full matrix. Maple identifies rank 3 exactly and the 
nontrivial eigenvalues in the interval [0.127,3.033]. 

Details for Example 14.41 The coefficients of this diagonally-implicit 4-stage 
method are 


/0.5544770574 0 

1.249724440 0.5544770574 

0.5816903621 -0.1016593761 
\ 1.591284531 0.2616450399 


/-0.1517881356 
0.423976329 
-0.2974468417 
\ 0.765740036 


-0.191898050 

0.56245777 

-0.421468111 

0.28293134 


0 

0 

0.5544770574 

-0.1745521833 

1.405289413 

0.163462193 

1.615228204 

-0.26943806 


0 

0 

0 

0.5544770574 

-0.061603226\ 
-0.14989630 
0.1036867507 
0.22076669 / 


\ 


where B has the eigenvalues = 1, 0.656, 0.397, 0.194. The difference matrix is 


/ -2.924587038 
0.2475622961 
-0.6719696591 
\ -0.391968580 


-2.026103023 

-0.163685467 

1.390690288 

-4.434411911 


4.236386571 

-1.071270925 

-0.3042527174 

1.433855269 


0.714303489 \ 
0.9873940939 
-0.4144679131 
3.392525220 / 


This method was found by a numerical search procedure in the Nordsieck formu¬ 
lation with a given slack Mu = diag(0, 0,1, 2) having the exact rank defect from 
Lemma 13.31 Transforming the computed weight matrices back to the original 
form gives 


/ 24.93687174 
36.6019187 
-41.1030151 
\-17.14663839 

/ 5.826768721 
17.59843495 
-15.32737123 
\ -9.04101962 


36.60191888 

104.2243208 

-93.94463413 

-50.66305601 

17.59843465 

58.3547424 

-47.2876058 

-30.63867973 


-41.10301520 

-93.94463411 

92.65849294 

44.55523963 

-15.32737112 

-47.28760605 

41.36332459 

24.30724431 


-17.14663842\ 
-50.66305601 
44.55523962 ’ 

24.93746522 / 

-9.0410195l\ 
-30.6386801 
24.3072443 ’ 

16.23315933 / 
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with smallest eigenvalues = 0.060 and = 0.038. With these matrices the test 
matrix M in (12.111) is full and has 2 numerical eigenvalues of small magnitude 
indicating rank 6 . However, the smallest eigenvalues is negative, —6 • 10“®. 
This confirms that without the transformation (|3.2I) semi-definiteness may not 
be verified reliably with real entries. 

The details of the reconstruction of this method from E by the Algorithm are 
as follows. Step 1 results in the matrix 


/ 1 
-0.9143733376 
1.226883712 
\ 0.1243293293 


0.9143733374 

1.557385735 

-1.704732078 

0.5657875861 


-0.8734437698 

-0.5451184012 

2.767561341 

-0.8827898483 


0.0868490964\ 
0.6055094088 
0.2225164874 
1.889057453 / 


Its triangular canonical form is H = UhEhU^^ with 


Lh 


Uh 


/ 1.803501131 
-2.878857141 

1.803501131 


\ 

1.236357981 - 

-0.9614900034 

1.803501131 


\ 0.048570432 

0.128230559 - 

-0.235314677 1 

.803501131; 

/I 0.1516606939 

-0.3187173423 

0.5591253567\ 


0.5468226337 

0.4769330471 

1.445134755 


V 

0.9431116674 

1.253027500 
2.559774554 ^ 

5 


and for the transformed seed matrix we obtain 


E' = Uh^EUh = 


/O 1 -1.411974971 
0 2 
0 


0.836862714\ 

0.226595572 

3 

0 y 


satisfying the condition 26^4 — 643(643 -I- 624 ) = 0 according to (14.41) . We note 
that known values which were distorted by rounding errors have been replaced 
(i.e., fill, the main diagonal of Lh = H' and the superdiagonal of E'). 

Details for Example 15.21 The coefficient matrices of this parallel 3-stage 
method are 



|. 

(K 

_ 1 _ 



1 ¥ 

¥ 


V Try 


\ 22 

22 

TT y 


and B has the eigenvalues = 1, 0.692 ± 0.172i. The original weight matrices are 


Z = H' ZH = 


145 

55 

24 

4 

841 

319 

80 

. 20 


484 

* 

15 


. /37 59 -91 \ 
w = e'^we =— * 137 -167 

\ * * 236 J 
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with minimal eigenvalues = 0.122 and = 0.021. Since the test matrix M in 
(I2.11|) is a full matrix with rational entries Maple is able to find rank 3 with 
nontrivial eigenvalues G [0.12,40.16]. 

Details for Example 15.31 For this parallel 4-stage method the coefficient ma¬ 
trices are 


/lO 

V 


16 


\ 

24 

30/ 


2 

25 

0 

52 

\ 



4^8 

3?^ 

21 

138 

35 

6 

W 

25 

35 

175 

116 

135 

165 

35 / 

35 

14 

14 


and B possesses the eigenvalues 1,3/5, ±1/5. The weight matrices are 


Z = H' ZH = 


w = e'^we = 


' 1299365 
63504 
* 

* 

* 

^ 743437 
79380 
* 

* 


329675 

2551§j^625 

32514048 

* 

* 


19656575 

6967296 

530265625 

8128512 

* 


221965 \ 
66 ®] 
^ m^5 

16j50?5 
190512 / 


3233059 

127008 

4662305 

63504 

* 


1101553 

42336 

93605 


1967785 

21168 


* 


530603N 


14112 

509231 

' MW 

6615 , 


with minimal eigenvalues = 0.109 and = 0.062. With these matrices the test 
matrix (|2.11l) can be computed. Its rank is 6 in accordance with the construction 
procedure and its nontrivial eigenvalues are contained in [0.021,194.1]. 
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